{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from ([this Jupyter notebook](https://nbviewer.jupyter.org/github/krischer/seismo_live/blob/master/notebooks/Computational%20Seismology/Wave%20Propagation%20%26%20Analytical%20Solutions/Greens_function_acoustic_1-3D.ipynb)) by Kristina Garina, Ashim Rijal and Heiner Igel ([@heinerigel](https://github.com/heinerigel)) which is a supplemenatry material to the book [Computational Seismology: A Practical Introduction](http://www.computational-seismology.org/),  additional modifications by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computation of Green's functions and seismograms for the acoustic wave equation\n",
    "\n",
    "In the previous lessons we derived the equations of motion to describe wave propagation in 3D, 2D and 1D elastic and acoustic media. Before solving the underlying partial differential equations numerically using finite-differences (FD), we should estimate some analytical solutions. \n",
    "\n",
    "This is not only useful to check if the FD codes contain any bugs, but also to get an idea of the accuracy of the numerical solution compared to the analytical."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1D Green's function  \n",
    "\n",
    "\n",
    "Let's start with a simple problem, like the equations of motion for a **1D acoustic medium** assuming a constant density model, which we derived in [Lesson 3](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/01_Analytical_solutions/3_Acoustic_medium.ipynb) and [4](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/01_Analytical_solutions/4_2D_1D_elastic_acoustic_approx.ipynb):\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 P}{\\partial t^2} - V_p^2\\frac{\\partial^2 P}{\\partial x^2} = f\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "If we introduce the **Dirac delta function**\n",
    "\n",
    "$$\n",
    "\\delta(x) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "\t\\infty &x=0 \\\\\n",
    "\t0 &x\\neq 0 \t\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "with the normalization condition \n",
    "\n",
    "$$\n",
    "\\int_{-\\infty}^{\\infty} \\delta(x)\\; dx = 1\n",
    "$$\n",
    "\n",
    "as a source term of the acoustic wave equation:\n",
    "\n",
    "\\begin{equation}\n",
    "f = \\delta(t-t_s) \\delta(x-x_s) \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "we call the solution of the 1D wave equation\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 G_1}{\\partial t^2} - V_p^2 \\frac{\\partial^2 G_1}{\\partial x^2} = \\delta(t-t_s) \\delta(x-x_s)\n",
    "\\end{equation}\n",
    "\n",
    "**Green's function** $\\mathbf{G_1(x,t;x_s,t_s)}$. This means that we place a source at $x = x_s$. The **source time function** describes the time-dependent behaviour of the source. In this case the source time function has an amplitude of $1\\; \\frac{Pa}{s^2}$ only at time $t = t_s$, otherwise the amplitude is zero. The pressure wavefield is recored at the receiver position x and time t.\n",
    "\n",
    "In the following derivation of the Green's function, we assume the special case of a source located at $x_s = 0\\; m$ and a source time $t_s = 0\\; s$, which simplifies eq. (1) to \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 G_1}{\\partial t^2} - V_p^2 \\frac{\\partial^2 G_1}{\\partial x^2} = \\delta(t) \\delta(x)\n",
    "\\end{equation}\n",
    "\n",
    "Furthermore, the P-wave velocity distribution in the sub-surface is constant:\n",
    "\n",
    "\\begin{equation}\n",
    "V_p(x) = V_{p0} = const. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "The Green's function can be calculated using different approaches. We try to find a solution in the Fourier domain. First, we apply a temporal **Fourier transform** \n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{f}(\\omega) = \\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty} f(t) e^{-i\\omega t} dt\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "where $\\omega$ denotes the angular frequency, to eq. (2):\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1}{2 \\pi} \\int_{-\\infty}^{\\infty}\\biggl\\{\\frac{\\partial^2 G_1(x,t)}{\\partial t^2} - V_{p0}^2\\frac{\\partial^2 G_1(x,t)}{\\partial x^2}\\biggl\\} e^{-i\\omega t} dt = \\frac{1}{2 \\pi} \\int_{-\\infty}^{\\infty}\\delta(t) \\delta(x) e^{-i\\omega t} dt \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Using the properties:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1}{2 \\pi} \\int_{-\\infty}^{\\infty}\\biggl\\{\\frac{\\partial^2 G_1}{\\partial t^2}\\biggl\\} e^{-i\\omega t} dt = -\\omega^2 \\hat{G}_1(x,\\omega) \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "and\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1}{2 \\pi} \\int_{-\\infty}^{\\infty}\\delta(t) e^{-i\\omega t} dt = \\frac{1}{2 \\pi} e^{-i\\omega 0} = \\frac{1}{2 \\pi}\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "we can get rid of the 2nd time-derivative on the LHS and the time-dependent $\\delta$ function on the RHS:\n",
    "\n",
    "\\begin{equation}\n",
    "-\\omega^2 \\hat{G}_1(x,\\omega) - V_{p0}^2\\frac{\\partial^2 \\hat{G}_1(x,\\omega)}{\\partial x^2} = \\frac{1}{2 \\pi} \\delta(x)\n",
    "\\end{equation}\n",
    "\n",
    "Next, we apply a spatial Fourier transform\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{f}(k) = \\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty} f(x) e^{-ikx} dx,\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "where $k$ denotes the wavenumber, to eq. (3):\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1}{2 \\pi} \\int_{-\\infty}^{\\infty}\\biggl\\{-\\omega^2\\hat{G}_1(x,\\omega) - V_{p0}^2\\frac{\\partial^2 \\hat{G}_1(x,\\omega)}{\\partial x^2}\\biggr\\} e^{-ikx} dx = \\frac{1}{4 \\pi^2}\\int_{-\\infty}^{\\infty}\\delta(x) e^{-ikx} dx\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "and using the same properties as for the temporal Fourier transform, we get:\n",
    "\n",
    "\\begin{equation}\n",
    "-\\omega^2 \\hat{G}_1(k,\\omega) + k^2 V_{p0}^2 \\hat{G}_1(k,\\omega) = \\frac{1}{4 \\pi^2}\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Solving for $\\hat{G}_1(k,\\omega)$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{G}_1(k,\\omega) = \\frac{1}{4 \\pi^2} \\frac{1}{V_{p0}^2 k^2 - \\omega^2} = \\frac{1}{4 \\pi^2 V_{p0}^2} \\frac{1}{k^2 - \\frac{\\omega^2}{V_{p0}^2}}\n",
    "\\end{equation}\n",
    "\n",
    "we have derived the **Green's function solution for the 1D acoustic wave equation in the frequency-wavenumber domain**.\n",
    "\n",
    "To get the time-domain solution, we first apply the inverse spatial Fourier transform to eq. (4)\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{G}_1(x,\\omega) = \\frac{1}{4 \\pi^2 V_{p0}^2}\\int_{-\\infty}^{\\infty} \\frac{e^{ikx}}{k^2 - \\frac{\\omega^2}{V_{p0}^2}}dk \\notag\n",
    "\\end{equation}\n",
    "\n",
    "This integral has two poles at $k = \\pm \\frac{\\omega}{V_p}$, so we have to integrate along a contour around the poles. For more details I refer to \n",
    "[S.W. Rienstra & A. Hirschberg (2017): An Introduction to Acoustics](https://www.win.tue.nl/~sjoerdr/papers/boek.pdf)\n",
    "\n",
    "The result is the solution in the **frequency-space domain**:\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{G}_1(x,\\omega) =  \\frac{e^{-i\\omega|x|/V_{p0}}}{4 \\pi i V_{p0} \\omega}. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "This will become useful to test the accuracy of our frequency domain finite-difference codes, which we will develop later in the lecture. For the transformation to time-domain we have to integrate around the pole at $\\omega=0$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{G}_1(x,t) =  \\frac{1}{4 \\pi i V_{p0}}\\int_{-\\infty}^{\\infty}\\frac{e^{i\\omega(t-|x|/V_{p0})}}{\\omega} d\\omega. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "This finally leads to the Green's function for the 1D homogeneous acoustic problem:\n",
    "\n",
    "\\begin{equation}\n",
    "G_1(x,t)=\\dfrac{1}{2V_{p0}}H\\biggl(t-\\dfrac{|x|}{V_{p0}}\\biggr),\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "where $H$ denotes the **Heaviside function**:\n",
    "\n",
    "$$\n",
    "H(x) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "\t0 &x<0 \\\\\n",
    "\t1 &x\\geq 0 \t\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "More generally, we can replace:\n",
    "\n",
    "\\begin{align}\n",
    "x &\\rightarrow x - x_s,\\nonumber\\\\\n",
    "t &\\rightarrow t - t_s,\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "and get:\n",
    "\n",
    "\\begin{equation}\n",
    "G_1(x,t)=\\dfrac{1}{2V_{p0}}H\\biggl((t-t_s)-\\dfrac{|x-x_s|}{V_{p0}}\\biggr),\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "So the 1D Green's function is a Heaviside function delayed by the traveltime between source and receiver. Note also that the absolute value of the offset $|x-x_s|$ implies that we have a wave propagating to the left and one propagating to the right. Let's plot the 1D Green's function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Define parameters\n",
    "        # velocity m/s \n",
    "        # distance from source\n",
    "        # length of seismogram (s)\n",
    "        # number of time samples\n",
    "        # time increment\n",
    "        # source time\n",
    "\n",
    "\n",
    "# Acquisition geometry\n",
    "            # coordinates of source\n",
    "\n",
    "\n",
    "\n",
    "            # coordinates of receiver\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Define time vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Calculating Green's function in 1D\n",
    "\n",
    "      # initialization G with zeros\n",
    "\n",
    "\n",
    "# Plotting Green's function in 1D\n",
    "plt.plot(time, G1)\n",
    "plt.title(\"Green's function for hom. 1D acoustic medium\" )\n",
    "plt.xlabel(\"Time, s\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D Green's function  \n",
    "\n",
    "Using the same approach as for the 1D medium, we can calculate the **2D Green's function** $\\mathbf{G_2(x,t;x_s,t_s)}$, which is governed by\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 G_2}{\\partial t^2} - V_{p0}^2 \\biggl(\\frac{\\partial^2 G_2}{\\partial x^2} + \\frac{\\partial^2 G_2}{\\partial z^2} \\biggl)= \\delta(t-t_s) \\delta(x-x_s) \\delta(z-z_s) \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "as \n",
    "\n",
    "\\begin{equation}\n",
    "G_2(x,z,t) = \\dfrac{1}{2\\pi V_{p0}^2}\\dfrac{H\\biggl((t-t_s)-\\dfrac{|r|}{V_{p0}}\\biggr)}{\\sqrt{(t-t_s)^2-\\dfrac{r^2}{V_{p0}^2}}} \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "with the source-receiver distance (offset)\n",
    "\n",
    "$r = \\sqrt{(x-x_s)^2+(z-z_s)^2}$\n",
    "\n",
    "Compared to the 1D Green's function, we have a damped Heaviside function due to the radiation characteristic of the infinite line source, introduced by the 2D approximation. Let's also plot the 2D Green's function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Calculation of 2D Green's function\n",
    "\n",
    "                    # initialization G with zeros\n",
    "\n",
    "           \n",
    "# Plotting Green's function in 2D\n",
    "plt.plot(time, G2)\n",
    "plt.title(\"Green's function for hom. 2D acoustic medium\" )\n",
    "plt.xlabel(\"Time, s\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "plt.xlim((0, tmax))\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3D Green's function  \n",
    "\n",
    "Finally, we can calculate the **3D Green's function** $\\mathbf{G_3(x,t;x_s,t_s)}$, which is governed by\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 G_3}{\\partial t^2} - V_{p0}^2 \\biggl(\\frac{\\partial^2 G_3}{\\partial x^2} + \\frac{\\partial^2 G_3}{\\partial y^2} +\\frac{\\partial^2 G_3}{\\partial z^2} \\biggl)= \\delta(t-t_s) \\delta(x-x_s) \\delta(y-y_s) \\delta(z-z_s) \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "as \n",
    "\n",
    "\\begin{equation}\n",
    "G_3(x,y,z,t) = \\dfrac{1}{4 \\pi V_{p0}^2 r}\\delta\\biggl((t-t_s)-\\frac{r}{V_{p0}}\\biggr) \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "with the source-receiver distance (offset)\n",
    "\n",
    "$r = \\sqrt{(x-x_s)^2+(y-y_s)^2+(z-z_s)^2}$\n",
    "\n",
    "So the 3D Green's function for the homogeneous acoustic medium is a Delta distribution delayed by the traveltime between source and receiver. For the computation of the 3D Green's function, we have to approximate the  $\\delta-$function. An example is the boxcar function\n",
    "\n",
    "$$\n",
    "\\delta_{bc}(x) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "\t1/dx &|x|\\leq dx/2 \\\\\n",
    "\t0 &\\text{elsewhere} \t\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "fulfilling the properties of the $\\delta$ function as  $dx \\rightarrow\\; 0$. This function is used to properly scale the source term to obtain correct absolute amplitudes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Calculation of 3D Green's function\n",
    "\n",
    "                                # initialization G with zeros\n",
    "\n",
    "                                # defining offset\n",
    "                                # defining amplitudes\n",
    "                                # time arrival\n",
    "\n",
    "# Plotting Green's function in 3D\n",
    "plt.plot(time, G3)\n",
    "plt.title(\"Green's function for hom. 3D acoustic medium\" )\n",
    "plt.xlabel(\"Time, s\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "plt.xlim((0, tmax))\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise\n",
    "\n",
    "Derive the Green's function solutions of the 2D elastic SH problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho\\frac{\\partial^2 u_y}{\\partial t^2} - \\frac{\\partial}{\\partial x} \\mu \\frac{\\partial u_y}{\\partial x} - \\frac{\\partial}{\\partial z} \\mu \\frac{\\partial u_y}{\\partial z} &= f_y\\nonumber \\\\\n",
    "\\end{align}\n",
    "\n",
    "and 1D elastic SH problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho\\frac{\\partial^2 u_y}{\\partial t^2} - \\frac{\\partial}{\\partial x} \\mu \\frac{\\partial u_y}{\\partial x} &= f_y.\\nonumber \\\\\n",
    "\\end{align}\n",
    "\n",
    "Assume a constant distribution of the shear modulus $\\mu(x,z) = \\mu_0 =  const. \\neq 0\\; Pa$ and density $\\rho(x,z) = \\rho_0 = const. \\neq 0\\; \\frac{kg}{m^3}$ in the subsurface.\n",
    "\n",
    "*Hint: To solve this problem, you do not have to apply any Fourier transform or integrate along contours around poles*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computation of seismograms\n",
    "\n",
    "In field data applications we can not excitate a source time function like a delta distribution, which would have a perfect white spectrum (all frequencies are excitated at once). As we will see later, some numerical problems arise when explicitly calculating Green's function with the time-domain finite-difference approach. \n",
    "\n",
    "Instead we have to rely on band-limited source signals. Seismograms for an arbritary source wavelet can be computed from the the Green's function. In the following example, the source wavelet consists of the first derivative of the Gaussian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Defining source time function\n",
    "           # Frequency (Hz)\n",
    "           # period\n",
    "           # defining t0     \n",
    "\n",
    "# Initialization of source-time function\n",
    "\n",
    "\n",
    "# Initialization of first derivative of gaussian\n",
    "\n",
    "\n",
    "# Plotting of source time function\n",
    "plt.plot(time, src)\n",
    "plt.title('Source time function')\n",
    "plt.xlabel('Time, s')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Excerise\n",
    "\n",
    "Compute seismograms $G_{seis}(x,t)$ for the 1D, 2D and 3D acoustic media, by a convolution of the Green's function $G(x,t;x_s,t_s)$ with the source wavelet $s(t)$:\n",
    "\n",
    "$$G_{seis}(x,t) = G(x,t;x_s,t_s) * s(t)$$\n",
    "\n",
    "Plot the resulting seismograms together with the Green's function solutions.\n",
    "\n",
    "*Hints:* \n",
    "* Use the NumPy function np.convolve. \n",
    "* How could you check if your implemented convolution is correct?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Computation of 1D seismogram\n",
    "\n",
    "# Convolution of Green's function with the 1st derivative of a Gaussian\n",
    "# COMPUTE YOUR SEISMOGRAM HERE!\n",
    "#G1_seis=\n",
    "\n",
    "# Plotting Green's function in 1D\n",
    "plt.plot(time, G1)\n",
    "plt.title(\"Green's function for hom. 1D acoustic medium\" )\n",
    "plt.xlabel(\"Time, s\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# Plotting convolved Green's function in 1D\n",
    "# PLOT YOUR SEISMOGRAM HERE!\n",
    "# plt.plot()\n",
    "plt.title('After convolution')\n",
    "plt.xlabel('Time, s')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.xlim (0, tmax)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Convolution of Green's function with the 1st derivative of a Gaussian\n",
    "# COMPUTE YOUR SEISMOGRAM HERE!\n",
    "#G2_seis=\n",
    "\n",
    "# Plotting Green's function in 2D\n",
    "plt.plot(time, G2)\n",
    "plt.title(\"Green's function in 2D\" )\n",
    "plt.xlabel(\"Time, s\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "plt.xlim((0, tmax))\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# Plotting convolved Green's function in 1D\n",
    "# PLOT YOUR SEISMOGRAM HERE!\n",
    "# plt.plot()\n",
    "plt.title('After convolution')\n",
    "plt.xlabel('Time, s')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.xlim((0, tmax))\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Convolution of Green's function with the 1st derivative of a Gaussian\n",
    "# COMPUTE YOUR SEISMOGRAM HERE!\n",
    "#G3_seis =\n",
    "\n",
    "# Plotting Green's function in 3D\n",
    "plt.plot(time, G3)\n",
    "plt.title(\"Green's function in 3D\" )\n",
    "plt.xlabel(\"Time, s\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "plt.xlim((0, tmax))\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "# Plotting convolved Green's function in 1D\n",
    "# PLOT YOUR SEISMOGRAM HERE!\n",
    "# plt.plot()\n",
    "plt.title('After convolution')\n",
    "plt.xlabel('Time, s')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.xlim (0, tmax)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We learned:\n",
    "\n",
    "* Calculation of analytical Green's function in the 1D, 2D, and 3D cases\n",
    "* How to define a source time function\n",
    "* Convolution of Green's function with source time function"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
